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Non-parametric strong lens inversion of SDSS J1004+4112 
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ABSTRACT 

In this article we study the weU-known strong lensing system SDSS J1004+4112. Not 
only does it host a large-separation lensed quasar with measured time-delay infor- 
mation, but several other lensed galaxies have been identified as well. A previously 
developed strong lens inversion procedure that is designed to handle a wide variety of 
constraints, is applied to this lensing system and compared to results reported in other 
works. Without the inclusion of a tentative central image of one of the galaxies as a 
constraint, we find that the model recovered by the other constraints indeed predicts 
an image at that location. An inversion which includes the central image provides 
tighter constraints on the shape of the central part of the mass map. The resulting 
model also predicts a central image of a second galaxy where indeed an object is visible 
in the available ACS images. We find masses of 2.5 x 10^'^ Mq and 6.1 x lO^'^ M© 
within a radius of 60 kpc and 110 kpc respectively, confirming the results from other 
authors. The resulting mass map is compatible with an elliptical generalization of a 
projected NFW profile, with — 58^23 arcsec and Cvir — 3.91 ± 0.74. The orientation 
of the elliptical NFW profile follows closely the orientation of the central cluster galaxy 
and the overall distribution of cluster members. 
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1 INTRODUCTION 

The gravitational deflection of light depends on both the 
luminous and dark matter present in the deflecting object, 
causing it to be an independent probe of the mass and possi- 
bly even the mass distribution of the deflector, i.e. the gravi- 
tational lens. Good alignment between a source, the lens and 
an observer, otherwise known as strong lensing, can cause 
several images of the source to be formed. Less perfect align- 
ment, or weak lensing, will not cause multiple images to ap- 
pear, but will still deform the image of the source somewhat. 
Multiple images and deformed images provide information 
about the mass distribution of the deflector and one can try 
to use these data to invert the lens, i.e. to determine its 
projected mass distribution. 

The lensing cluster SDSS J1004-I-4112 was revealed by 
the presenc e of a multiply imaged quasar as reported by 
llnada et all (|2003l ). The lensing system was first identified 
as a quadruply imaged quasar, but later a fifth cen tral image 
of the quasar was detected by 'inada et al.' ('2005') and spec- 
troscopically confirmed by[lnada ot al. (2008). Three multi- 
ply imaged galaxies were identified in HST/ACS images by 



ISharon et all l|2005l ) and time delay info rmation for three 
of the quasar images was measured by iFohlmeister et al] 
(|2008t ). improving the earlier rep orted time delay betwee n 
the two closest quasar images in IFohlmeister et al.l (|2007h . 
This work did not only only invalid ate earlier p r opose d 
models of the l e nsing system (e.g. lOguri et al.l l|2004h . 
IWiUiams fc Sahal i20o4 )). which predicted shorter time de- 
lays, it also confirmed that microlensing is the cause of 
the strange magnificatio n patterns in the qu asar images, 
resent both in op tical l|Richards et al.l 120041 ) and X-ray 
Lamer et al.l |2006| ) measurements. With its separation of 
14 arcsec, the multiply imaged quasar in SDSS J1004-I-4112 
has held the record for being the widest lensed quasar for 
a number of years. The discovery of SDSS J1029-I-2623, a 
multiply imaged q uasar with a separation of over 22 arcsec 
(|lnada et al.l |2006| ) broke this record recently. The statis- 
tics of multiply imag ed quasars by clusters are studied in 



tics 01 multiply imag e 
iHennawi et all ||20071 ). 
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In a strong lensing scenario, various kinds of informa- 
tion can be available, all encoding some information about 
the projected mass distribution of the lens. Not only does 
one have positional information of images of the same source, 
but it is also possible that magnification information or time 
delay information is present. Even the absence of images 
in certain locations can provide constraints on the mas s 
distribution. In previous works l|Liesenborgs et al.l l|2006h . 
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iLiesenborgs et alj l|2007l ') and lLiesenborgs et all (|2008bl )') we 
described a flexible, non-parametric method for strong lens 
inversion. In this article, we shall apply this procedure to 
SDSS J1004+4112 and compare our results with other find- 
ings about this system. 

In section [2] we will review the non-parametric inversion 
method described in previous works, and discuss some mod- 
ifications and extensions. We shall apply this method to the 
gravitational lensing system SDSS J1004-I-4112 in section|3l 
Finally, the results of this inversion will be discussed in sec- 
tion U Unless noted otherwise, uncertainties mentioned in 
this article specify a 68% confidence level. 



2 INVERSION METHOD 

2.1 Lensing basics 

Image formation in gravitational lensing is most often de- 
scribed by the lens equation, which relates points G on the 
image plane, describing what can be seen because of the de- 
flection of light, to points /3 on the source plane, describing 
what one would see without the lens effect: 

OiG) ^ e - ^&{e). (1) 

The angular diameter distances D^s and Ds measure the 
distance between lens and source and observer and source 
respectively, and depend on the redshifts of source and lens. 
The actual bending of light rays is stored in a , the deflection 
angle. 

The light travel time from source to observer is de- 
scribed by the time delay function: 

m /3) = i±il^ Q(0 - /3)^ - ^ . (2) 

Here, Zd describes the redshift of the gravitational lens, with 
corresponding angular diameter distance Dd- The gradient 
of the lensing potential tp is related to the deflection angle a 
in such a way that the stationary solution of the time-delay 
function, i.e. Vet = 0, again yields the lens equation. 

If the gravitational lens effect causes a single source to 
be seen as several distinct images, different light travel times 
will give rise to a time delay. Denoting /3 the source position 
and 9i and 02 two corresponding image positions, the time 
delay between the two images is given by 

Ati2 = i(6li,/3)-t(6»2,/3). (3) 

If the source brightness is time- variable, similar brightness 
variations will be seen in the images at different times, and 
this time delay may be measured. 

For more detailed information about the gravitational 
lensing formali s m, th e interested reader is referred to 
ISchneider erall l| 19921 ). 

2.2 Genetic algorithm based inversion 

As described in ILiesenborgs et aL the inversion 

method we propose requires the user to specify a square- 
shaped region in which the procedure should try to recon- 
struct the projected mass density E of the lens. At first, this 
region is subdivided into a number of smaller squares in a 
uniform way, and to each square a projected Plummer sphere 



(|Plummej|l91ll ) is assigned. A genetic algorithm then looks 
for appropriate weights of these basis functions to construct 
a first approximation of the mass distribution. Using this 
first solution, a new grid is created in which regions con- 
taining more mass are subdivided further and the genetic 
algorithm again tries to determine appropriate weights for 
the associated basis functions. This iterative scheme can be 
repeated until the added resolution no longer considerably 
improves the fit to the data. 

It is clear that the genetic algorithm mentioned above 
is the core of the inversion procedure. A genetic algorithm is 
an optimization strategy inspired by the Darwinian theory 
of evolution. An initial population of random trial solutions 
is evolved into solutions which are better adapted to the 
problem under study. To create each new generation, trial 
solutions are combined, cloned and mutated, and in doing 
so, selection pressure must be applied: trial solutions which 
are considered to be more fit, should create more offspring. 
Not only is it possible this way to create solutions which are 
optimized with respect to a single criterion, but so-called 
multi-objective genetic algorithms allow several fitness mea- 
sures to be optimized at the same time. A detailed account 
of genetic algorithms and multi-objective genetic algorithms 
in particular, can be found in Deb (200i|). The different fit- 
ness criteria that we use will be discussed below. 

The origin al procedur e as describ ed in 

ILiesenborgs et all l|2006h and ILiesenborgs et all (|2007h 
has a shortcoming which is illustrated in Fig. [T] The left 
panel of this figure shows a mass map which consists of 
relatively small density peaks on top of a sheet of constant 
density. When our procedure is used to reconstruct the 
projected density of this lens, it fails in a quite dramatical 
manner (center panel). There are two causes of this unde- 
sirable behavior. First, the algorithm will have to try to 
mimic the effect of a mass sheet using the Plummer basis 
functions which is a rather difficult task, depending on the 
amount of constraints available. The second problem is that 
the subdivision scheme will be less effective. Since the mass 
sheet holds most of the mass, the subdivision procedure will 
not be successful at refining the grid in the central region. 

To solve these problems, we have made it possible for 
the user to specify that the algorithm should also search for 
a mass sheet. In effect, we have added a mass sheet as a basis 
function for which the genetic algorithm needs to determine 
the weight. The subdivision scheme can then inspect the 
mass density relative to this sheet of mass so that the regions 
of interest can again be reconstructed with a finer resolution. 
The right panel of Fig. [T] illustrates how much this can help 
to improve the reconstruction. 

The entire inversion procedure can be repeated a num- 
ber of times to create a number of solutions which are com- 
patible with the input constraints. Using such a set of so- 
lutions one can inspect the average, which highlights the 
common features of the mass maps, and one can calculate 
the standard deviation, revealing the areas in which solu- 
tions tend to disagree about the exact shape of the projected 
density. 
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Figure 1. Left panel: true projected mass density of a lens used to test the inversion procedure. The mass distribution consists of a few 
relatively small perturbations on top of a sheet of mass. Center panel: when the original procedure is applied to the images produced by 
the input lens, it is not successful in creating an acceptable mass map (see text). Right panel; when a sheet of mass is added as a basis 
function, the algorithm again is able to create acceptable reconstructions of the projected density. 



2.3 Fitness criteria 

2.3.1 Positions 

In strong lens inversion, an obvious set of constraints is the 
information about multiply imaged sources. If the true mass 
distribution of the gravitational lens were known, using it to 
project the images of a single source back onto the source 
plane, would result in a single consistent source shape. If an 
incorrect lens is used, the lens equation will project each im- 
age back onto different regions in the source plane. The first 
fitness measure is therefore the amount of overlap between 
the back-projected images of each source. 

Each back-projected image of a single source is sur- 
rounded by a rectangle and the distances between corre- 
sponding corners of the rectangles are used to calculate a 
fitness value. If corresponding points in the images can be 
identified, they too can be included in the fitness measure. 
Note that in calculating such distances, the estimated source 
size is used as the length scale. This a voids over-focussing 
the images fsee lLiesenborgs et al.l (|2006l )') which is even more 
important when a mass sheet is included as a basis function: 
solutions with a considerable mass sheet will automatically 
project the images onto a smaller region in the source plane. 



2.3.2 Null space 

Using only the first criterion, the genetic algorithm evolves 
towards solutions for which the back-projected images over- 
lap. However, it is also possible that other regions of the im- 
age plane are projected onto the same region in the source 
plane. If this is the case, the suggested solution would pre- 
dict additional images. In situations where there are clearly 
no other images present, one would like to use this so-called 
null space as an additional constraint. 

To do so, the null space is subdivided into a number 
of triangles, and the trial solution under study is used to 
project these triangles onto the source plane. Then, the 
amount of overlap between each triangle and the current 
estimate of the source shape is calculated and used to con- 
struct a null space fitness measure. The envelope of the back- 
projected images is used to estimate the source shape. More 



detailed information about the u se of the null-space can be 
found in iLiesenborgs et al.l (|2007l ). 



2.3.3 Critical lines 

In many cases it is obvious that images are not intersected 
by a critical li ne, i.e. that all po i nts of an image have the 
same parity. In ILiesenborgs et al.l (I2OO8IJ ) we described how 
this information was used to avoid the genetic algorithm 
being trapped in a sub-optimal region of the solution space, 
where an image does get intersected by a critical line. The 
solution that was used simply calculated the sign of the 
magnification at several points inside an image, and this 
was used to construct a fitness measure which penalizes im- 
ages in which the sign changes. While this worked well in 
the case of CL 0024-1-1654, applying the same method to 
SDSS J1004-h4112 was far less successful. 

Fig. [5] illustrates the problem. In the left panel, the 
black regions mark two images of a single galaxy, and the 
points in each image should all have the same parity. For the 
constructed solution, the critical lines are shown and they 
clearly do not intersect the input images, meaning that no 
parity changes will be present in an image and that the so- 
lution will not be penalized. When the proposed mass map 
is used to project the images back onto the source plane, the 
situation in the right panel arises. Clearly, when the enve- 
lope of the back-projected images (grey area) is considered, a 
caustic does intersect this region and correspondingly when 
this shape is used to predict the images, a critical line will 
intersect an image as can be seen in the left panel. By not 
specifying precisely what type of solution one is interested 
in, the existing criterion can easily lead the genetic algorithm 
towards a sub-optimal reconstruction. 

Instead of calculating the magnification information at 
the location of the images, the value of the magnification is 
now calculated on a relatively coarse grid covering the re- 
gion of interest. This is used to create a rough estimate of 
the critical lines, which in turn are projected onto the source 
plane to provide an estimate of the caustics. The intersec- 
tion of the caustics with the source shape is calculated and 
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Figure 2. Illustration of the problem with the original fitness measure to penalize situations in which a critical line crosses an image. 
Suppose two input images (left panel, black) are known not to be intersected by a critical line. The critical lines of a certain trial solution 
indeed do not intersect the input images, so all the points in the input images will have the same parity. However, when the images are 
projected onto the source plane (right panel), the envelope of both images is in fact intersected by a caustic, causing a critical line to 
intersect the current prediction of the images (left panel, grey). 



the total length is used as a fitness measure, a lower value 
indicating a better fitness. 



2.3.4 Time delay information 

When time delay information is available for a number of 
images of a single source, one would like to use this infor- 
mation to constrain the allowed region in the solution space 
even further. By calculating the lensing potential at the im- 
age points for which time delay information is available, in 
principle equation ((2)| can be used to compare the predicted 
time delays with the observed ones. However, to do so, one 
needs to know the position (3 of the source. While the source 
position may be estimated once a good overlap of the images 
has been reached, this is in general not possible while the 
genetic algorithm is still evolving, and certainly not near the 
start, when the trial mass maps are still quite random and 
the images are projected onto very different regions. 

Having tested a number of possible fitness measures, 
we found that the following one works very well. Suppose 
that there are images Gi with corresponding points in the 
source plane /3^. It is possible that time delay information is 
not available for all images, so let us call T the set of image 
indices for which time delay information is at hand. The 
measured time delay between image i and j will be called 
Atobs.ij- The fitness measure is then given by: 



EEEE 

ieT j£T fe=l i = l 



[t(0i,/3fc) -f(6>,,/3,)] - Atobs,», 



At, 



ohs, ij 



(4) 



Again, a lower value implies a better fitness of the trial so- 
lution. 



3 APPLICATION TO SDSS J1004+4112 
3.1 Multiple image systems 

Fig. [3] shows the image systems that were used in the in- 
version of SDSS J1 004-I-4112, using the same labeling as 
ISharon et all (120051 ). There are five spectroscopically con- 
firmed images of a quasar at redshift 1.734, labeled Ql- 
Q5. Corresponding to the time delay measurements of 
iFohlmeister et all (|2008l ). we used a time delay of 40.6 days 
between Q2 and Ql, and a time delay of 821.6 days between 
Q3 and Ql. No magnification information was used, as the 
quasar image magnifications are influenced by microlensing, 
introducing a large uncertainty. The po sitions of the quasa r 
images were set to those reported in llnada et al. I l|2005l ). 
Four, possibly five images are present of a galaxy at redshift 
3.332, l abeled A1-A5 , with image A5 being marked as uncer- 
tain bv lSharon et al.l (|2005l ). The third system used consists 
of two images of a galaxy at redshift 2.74, marked B1-B2. 
Note that another galaxy with two images was identified in 
the aforementioned work, but because of its unknown red- 
shift, it was not used in the inversion. Angular diameter 
distances were calculated in a flat cosmological model with 
Ho = 71 km s"^ Mpc"^ D.^ = 0.27 and = 0.73. Us- 
ing the redshift information described above, this fixes the 
Dds/Ds ratios for the lensing systems, which is required in- 
put information in our method. 



3.2 First inversion 

Since image A5 was marked as uncertain, the first inversion 
does not include it. The algorithm was instructed to look for 
mass in a square region, 35 arcsec wide, roughly centered on 
image Q4. The null space fitness measure was based on a 
square region, 60 arcsec wide, subdivided into a 64 by 64 
grid. For each source, the image regions were excluded from 



Non-parametric inversion of SDSS J1004+4112 5 






X (Arcsec) 



a 
< 





X (Arcsec) 



Figure 4. Left panel: when the input images are projected back on to the source plane using the average of 28 individual solutions, 
these source positions are obtained. Galaxy A is surrounded by a dashed rectangle, galaxy B by a dotted one. The caustics correspond 
to the redshift of the quasar. Right panel: when the sources and caustics of the left panel are used to predict the images and critical lines 
using the average solution, this configuration is obtained. 
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Figure 5. Left panel: average mass map of 28 individual solutions when image A5 is not taken into account. Most of the mass is found 
to coincide with the region of the BCG. The critical density was calculated at the redshift of the quasar. Right panel: standard deviation 
of the individual solutions, showing that the precise distribution near the center of the cluster is somewhat uncertain. 



the null space, and for systems A and B, the central clus- 
ter region was excluded as well, allowing the algorithm to 
predict the locations of the central images of these systems. 
The null space is a relatively large region, but this avoids 
the introduction of unnecessary substructure at the edge of 
the mass reconstruction region, that would cause images to 
appear at larger distances. The critical line fitness was based 
on a square shaped region, 40 arcsec wide, subdivided into 
a 64 by 64 grid. After eac h inversion, a final izing step was 
performed, as described in iLiesenborgs et"aL 1 |2008b ). This 



causes some minor modifications to be made to the mass 
map, to improve the positional and time delay fitness mea- 
sures. In the same work we described how mass could be 
redistributed without affecting any of the observable prop- 
erties and demonstrated this on the obtained mass map for 
CI 0024-1-1654. In this work however, no explicit mass redis- 
tribution step is performed. 

The average solution of 28 individual inversions predicts 
the source positions and caustics shown in the left panel 
of Fig. |4] The source position of galaxy A is marked by a 
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Figure 3. The multiple image systems which are used in 
the inversion of SP SS J1004+4112, using the same labeling as 
ISharon et alj ||2005| ) (north is up, east is left). Five images of a 
quasar (Q1-Q5) are available, as well as four, possibly five images 
of a galaxy marked A1-A5, and two images of a second galaxy 
marked B1-B2. Between Bl and Q3 and to the left of B2 are tw o 
images of a third galaxy marked C1-C2 in lSharon et al 
but this system was not used as no redshift is currently available. 



dashed rectangle, the position of galaxy B is marked by a 
dotted one. When these sources and the reconstructed lens 
are used to predict the image configurations, the result in 
the right panel of the same figure is obtained. The critical 
lines and caustics in these figures are calculated for the red- 
shift of the quasar. The mass map itself is shown in the left 
panel of Fig. (5] with most of the mass in the same region 
as the brightest cluster galaxy (BCG). The standard devia- 
tion of the individual reconstruction can be seen in the right 
panel of the same figure, showing that the precise distribu- 
tion of mass in the central region differs between the indi- 
vidual reconstructions. Fig. [6] shows the average profile and 
its standard deviation. The large core clearly differs from 
the NFW-like behavior that one might expect. 

When inspecting the right panel of Fig. U one sees that 
the average solution predicts central images of galaxies A 
and B. The predicted position of the central image of galaxy 
A coincides with the location of image A5, although the pre- 
dicted shape is far less extended. Fig. [7] shows the central 
region of the cluster, after subtractin g the central clu ster 
members using the GALFIT software (|Peng et aLlbOOd ). In 
each of the filters, one can clearly see the central image of 
the quasar in the upper-left region. Image A5 can clearly be 
seen in the F555W and F814W images. Since the other con- 
straints predict a central image of galaxy A at this location 
and since it indeed resembles a mirror image of Al, we feel 
confident that this is in fact the central image of galaxy A. 
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Figure 6. The circularly averaged profile of the inversions when 
image A5 is disregarded, together with the standard deviation. 



Figure 7. The central part of the cluster after removing the 
contribution of the central cluster members using GALFIT. The 
central quasar image can clearly be seen in each filter, in the 
upper left part of the image. Below and to the left of it, image 
A5 can be seen in the F555W and F814W images. More to the 
right, an extra object can be seen, where the inversion predicts 
the central images of galaxies B and C. 
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Figure 8. When using the model resulting from the second inver- 
sion to project the galaxy images back onto their source planes, 
these images are obtained. Note that image A4 is not shown here, 
as it is occluded by a cluster galaxy. The size of galaxy A is ap- 
proximately 4 kpc, the size of galaxy B is approximately 2.5 kpc. 
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Figure 9. Left panel: when the average of 28 individual solutions is used to reconstruct the source plane when image A5 is included 
as a constraint, this result is obtained. The dashed box again indicates galaxy A, the dotted one galaxy B. Right panel: the sources 
and caustics in the left panel correspond to these images and critical lines. In this case, the central image of galaxy A is indeed more 
elongated. The critical lines and caustics again correspond to the redshift of the quasar. 
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Figure 10. Left panel: average mass density of the 28 individual solutions. When image A5 is included, the central region clearly needs 
to be much steeper. Right panel: standard deviation of the individual solutions. The precise mass distribution in the central region differs 
somewhat among the reconstructions. The critical density again corresponds to the critical density at the redshift of the quasar. 



3.3 Second inversion 

Including the central image of galaxy A will provide ad- 
ditional information that will lead to a different inversion 
since its true shape is different from the one predicted by 
the first inversion. For this reason, a second inversion was 
performed in which image A5 was added as an observational 
constraint. The rest of the constraints are the same as in the 
first inversion. Fig. [9] shows the source and image configura- 
tions obtained in this case, using the average solution of 28 
individual reconstructions. The central image of galaxy A is 



now clearly more extended than in the first inversion. When 
the images of galaxies A and B are projected back onto their 
source planes, the source shapes in Fig. |8] are reconstructed. 
The back-projected images of each source clearly resemble 
each other, illustrating that a good positional fitness has 
been achieved. The estimated size of galaxy A is approxi- 
mately 4 kpc, the size of galaxy B approximately 2.5 kpc. 

The effect of the inclusion of image A5 can best be 
seen in the average mass map, as shown in the left panel of 
Fig. 1101 Now, the mass distribution has clearly become much 
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Figure 12. Left panel: average profile and standard deviation of the resulting mass distributions. The dashed line shows the best fit 
NFW profile. Right panel: when only the mass density at the location of the images is taken into account, this is the resulting best fit 
NFW. The center of the profile lies very close to Q5, as does the center of the B CG. The orientati on is very similar to that of the BCG 
(dashed line), and corresponds to the general alignment of the cluster members ijOguri et al.ll20o3l . 




Figure 11. The average solution resulting from the second inver- 
sion is shown as a contour map on top of the ACS image. Most 
of the mass clearly lies in the same area as the central cluster 
members. The mass peak in the north-west part of the figure is 
not significant, as it can easily be redistributed. The dashed line 
indicates the orientation of the BCG. 

steeper in the central region, although some disagreement 
still remains between the individual solutions (right panel). 
A comparison with the visible matter can be seen in Fig. 1111 
The effect on the mass density can also be clearly seen in 
the circularly averaged profile, shown in the left panel of 
Fig. 1121 It would definitely be interesting to see how much 
the resulting mass map resembles a NFW distribution. 



The NFW density profile (|Navarro et all 119961 ) is de- 
scribed by: 

PNPw(r) = , , ../'^ — j-^, (5) 
(r/rs)(l + r/rs)2 

in which pa is a density scale factor and Vs is a characteristic 
radius. The density scale can be expressed in terms of Cvir, 
which relates to the virial radius rvir through rvir = CvirJ's- 
The virial radius itself is defined as the radius within which 
the mean density equals Avir times the mean matter density 
at the redshift of the halo. This virial overdensity Avir stems 
from the spherical collapse model, and fo r a fiat cosmologi- 
cal model it can be approxim ated by (e.g. iBrvan fc Normar] 
l|l998t ). lBullock et al.1 (120011 )) 

IStt^ + 82x - 39x^ 

' 

in which x — Q{z) — 1 and ^l{z) is defined as the ratio of the 
mean matter density to the critical density. Through lens 
inversion one recovers the projected density: 

pNFw{R,z)dz, (7) 

- oo 

for which an analy t ical e xpression can be calculated (e.g. 
IWright fc BrainerdI (|2000l )). 

Naively performing a fit of the profile in the left panel 
of Fig. [12] to a projected NFW profile, yields the best fit 
profile described by the dashed line in the same figure. One 
then finds Ts = 41.2^i['3 arcsec, and Cvir ~ 5.37^q j2- 
th ough this s e ems t o correspond well to the values found 
by lota et al.l (j2006l ). who reported rg = 39lg'^ arcsec and 
Cvir = 6.115^ 2 (90% confidence) based on Chandra X-ray 
observations, the uncerta inties found in th i s way are far 
too low. As explained in iLiesenborgs et al.l (|2008bl ). using 
the monopole degeneracy it is possible to redistribute the 
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Prediction 




CASTLES^ 




12005 


F2008 


Image 




F160W 


F555W ~ 


F814W 






Qi 


1 


1 


1 


1 


1 


1 


Q2 


1.03 ±0.38 


0.6486 


1.0864 


1.3428 


0.732 


0.724 


Q3 


0.54 ±0.19 


0.4487 


0.4529 


0.4656 


0.346 


0.592 


Q4 


0.29 ±0.11 


0.3191 


0.6138 


0.2489 


0.207 




Q5 


0.032 ±0.029 


0.0114 


0.00024 


0.0047 


0.003 





Table 1. The predicted flux ratios of the quasar images, compared to data from the CASTLES project, llnada et al ] l|2005h and 
iFohlmeister et al.l | |2008l ) respectively. Note that only in this last work, the combined effect of the intrinsic variability of the source 
and the time delay has been taken into account. The general trend of the predicted values matches the observations, even though no 
magnification information was used in the inversion. The uncertainties show that this non-parametric inversion method can create a wide 
variety of flux ratios, even without having to consider microlensing. 



mass in between the images, without affecting any of the ob- 
servable properties of the lensing system. This means that 
the uncertainty of the circularly averaged profile is actu- 
ally much larger than obtained by simply calculating the 
standard deviation of the individual profiles. In turn, this 
translates to larger uncertainties on the parameters of the 
fit. 

Since the mass distribution in between the images is 
not well constrained, it is interesting to see how much the 
density at the location of the images themselves constrains 
the NFW parameters. First, we calculated the average den- 
sity and its standard deviation at the location of each image. 
Then, an elliptical generalization of Enfw was fitted to these 
data points. An axis ratio / was introduced in the projected 
NFW profile by setting R = {fx^+y'^/ff'^ in equation Q. 
We prefer this substitution over R= {x^ + (y/qfy^^ that 
would correspond to an axisymmetric NFW instead of a tri- 
axial one, because the circularly averaged profile in the first 
case corresponds closely to the profile of a symmetric NFW 
with the same rs and Cvir parameters. This allows the ob- 
tained values to be compared directly to fits to the circularly 
averaged profile. After fitting the elliptical generalization of 
Snfw, the values rg — 58^13 arcsec and Cvir ~ 3.91 ± 0.74 
are obtained. The best fit NFW is shown in the right panel 
of Fig. 1121 Its orientation corresponds to that of the BCG 
and to the general configurati on of the cluster members as 
reported in loguri et"al] (|2004l '). 

When calculating the total mass within 60 kpc, corre- 
sponding to the region of the quasar images, and 110 kpc, 
the region bounded by the images of galaxy A, we find re- 
sults of 2.5 X 10^^ Mq and 6.1 x 10^^ Mq r espectively. These 
values can be compared to the findings of I Williams fc Saha 

2004) . who also find 2.5 x 10^^ Mq, and of ISharon et al. 

20051) , who find 6 x 10^^ Mq- This illustrates nee more 



that th e mass within the im ages is well constrained. 

In ISharon et"ai] l|2005l ) a lens model was used to pre- 
dict the redshift of galaxy C, of which the two images lie 
between Bl and Q3, and to the left of B2 respectively (see 
Fig. [Sjl . Doing the same using the average model discussed 
above, we find that the back-projected images nearly over- 
lap for a Dds/Ds ratio of 0.64, corresponding to a redshift of 
3.35, slightly higher than the reported redshift of 2.94. Af- 
ter the inversions were completed, we have learned that the 
authors of the aforementioned work have now spectroscopi- 



http: / /www. cfa.harvard.edu/castles/ 



cally confirmed the redshift of galaxy C to be 3.288 (private 
communication) . 

The right panel of Fig. |9] contains a prediction for the 
central image of galaxy B, lying to the right of image A5. 
Inspecting Fig. [7| again, there indeed seems to be an ob- 
ject at that location, which is especially clear in the F435W 
and F555W filters. It is important to note however that the 
model also predicts that the central image of galaxy C men- 
tioned above, is located at almost the same location as the 
central image of galaxy B. For this reason, the object that 
can be seen in Fig. [T] is possibly a superposition of the cen- 
tral images of these two galaxies. 

The predicted flux ratios for the quasar system - rel- 
ative to the flux of Ql - are shown in table [T] and are 
compared to the flux ratios from other works. Although 
no magnification information was used in the inversion, the 
general trend of the predictions matches the observations. 
Also note that the relatively large uncertainties show that 
the non-parametric technique can accomodate a wide num- 
ber of flux ratios, without taking microlensing into account. 
Finally, the model presented here predicts a time delay of 
slightly over 1300 days between images Ql and Q4 of the 
qu asar. This is still c o nsiste nt with the constraint presented 
in IFohlmeister et ahl (|2008h which specifles that this delay 
should be over 1250 days. The Q1-Q5 time delay is predicted 
to be of the order of 1900 days. 



4 DISCUSSION AND CONCLUSIONS 

In this article we have applied a previously developed strong 
lens inversion method to the case of SDSS J1004+4112. The 
constraints used include time delay information, positional 
information and null-space information, all handled well us- 
ing a multi-objective genetic algorithm. 

The system under study only provides a few sources 
at different redshifts, which, in principle, still allows a gen- 
eralized version of the m ass sheet or steepness degeneracy 
(|Liesenborgs et al]|2008al ) . It is for this reason that the avail- 
able time delay information is of particular importance here, 
as it directly breaks the degeneracy. The fact that the degen- 
eracy is broken well can be seen in the low dispersion in the 
outer regions of the surface density (right panels of Figs. [5] 
and 1101) which is of the order of E/Ecr ~ 0.05, indicating 
that in our extended version of the genetic algorithm a sim- 
ilar mass sheet basis function is found in each individual 
reconstruction. It is interesting to compare the mass map of 
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the se cond inversion to the mass map obtained bv lSaha et al.l 
l|2007l ). The outer contours of their reconstruction show a re- 
markably circular structure, causing a similar eflFect as the 
mass sheet basis function used in our work. The contour 
steps in that figure would correspond to E/Ecr ~ 0.22, in- 
dicating that a similar mass density will be found near the 
edges of image system A as in our work. 

Note that in the reconstruction of the projected mass 
density, relatively large structures seem to exist to the north 
and south of images A3 and A4. As already suggested by the 
large associated standard deviations, one should not place 
much confidence in the displayed shape of these features, as 
the mass in those regions can easily be redistributed with- 
out affecting any of the observable prope rties of the lensing 
system using the monopole degeneracy (jLiesenborgs et aP 
l2008bl ). For the same reason it is extremely difficult to make 
reliable statements about the nature of substructure that 
may be present near the cluster center. One can only hope 
to make reliable predictions about the projected density at 
the location of the images themselves, illustrating the need 
for lenses with many multiply-imaged systems. Furthermore, 
to probe the core regions of clusters, central images are of 
particular importance as is nicely illustrated by the differ- 
ence in profiles between the two inversions shown in this 
article. 

When studying the constraints provided by the density 
at the image locations, we find that the resulting best fit 
NFW bears great resemblence to th e general clus t er con - 
figuration. As is often the case (e.g. iKeeton et al.l (|l998l )) 
the fit has a very similar orientation as that of the central 
galaxy, which in this case also follows the gene ral distribu- 
tion o f the cluster galaxies. In a recent study, lOguri et al.l 
l|2009l ) discussed the fact that lensing clusters are often over- 
concentrated. Although the circularly averaged profile in- 
deed suggests that this may be the case in this cluster as 
well, the more reliable two-dimensional fit yields an estimate 
of the concentration which is compatible with the expected 
value Cvir ~ 4. 

The method described and applied in this article is a 
non-parametric one, in the sense that no predefined shape 
for the matter distribution is used to fit the data. This is 
done by arranging a large number o f Plummer basi s func - 
tions on a grid. In a recent article, iJuUo fc KneTbl (|2009l ) 
made the interesting point that when basis functions over- 
lap, the introduced correlation reduces the effective number 
of degrees of freedom, making such a non-parametric inver- 
sion less underconstrained than it appears at a first glance. 
In any case, non-parametric methods can certainly help to 
explore a larger portion of the solution space, helping one to 
obtain a less biased look at the possible mass distributions. 
As with any method, one must be cautious about interpret- 
ing the results, since degeneracies can greatly enhance the 
uncertainties involved. 
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